1 Libraries

  library(package = "tidyverse")
## ── Attaching packages ───────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.1.0       ✔ purrr   0.3.2  
## ✔ tibble  2.1.1       ✔ dplyr   0.8.0.1
## ✔ tidyr   0.8.3       ✔ stringr 1.4.0  
## ✔ readr   1.3.1       ✔ forcats 0.4.0
## ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
  library(package = "lubridate")
## 
## Attaching package: 'lubridate'
## The following object is masked from 'package:base':
## 
##     date
  library(package = "rlist")
  
  library(package = "ncdf4")
  library(package = "ncdf4.helpers")
  
  library(package = "PCICt")
  
  library(package = "reshape2")
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths

2 Input Data File

# directory/URL root

  directory = "/projects/BIG_WEATHER/GENS_ERROR_RT/triangle_archives/"
  directory = "http://kyrill.ias.sdsmt.edu:8080/thredds/fileServer/BWW_GENS/CI_STAT/"

  time_range = "2015-06-01_00_to_2019-03-31_00"
  
  region = "WRFRAP"
  
  RdataFile = str_c(directory,
                  "gens_03_ensemble__",
                  "T2M_MSLP_M10_SPCH2M_ISOHGT_U10_V10_FRICV_GUST",
                  "__error__",
                  region,
                  "__",
                  time_range,
                  ".RData",
                  sep = "")

  load(file = url(RdataFile)) 
  
  Time     = unique(CI_Ensemble_Stats$Time)
  Fx_Hour  = unique(CI_Ensemble_Stats$Fx_Hour)
  Variable = unique(CI_Ensemble_Stats$Variable)
  Height   = unique(CI_Ensemble_Stats$Height)
  
  CI_Ensemble_Stats$Month   = month(CI_Ensemble_Stats$Time)

  CI_Ensemble_Stats = CI_Ensemble_Stats %>% 
    mutate(Quarter = case_when(((Month == 12) | (Month == 01) | (Month == 02))  
                               ~ "DJF",
                               ((Month == 03) | (Month == 04) | (Month == 05))
                               ~ "MAM",
                               ((Month == 06) | (Month == 07) | (Month == 08))
                               ~ "JJA",
                               ((Month == 09) | (Month == 10) | (Month == 11))
                               ~ "SON") )

3 Triangle Plots

3a 2-m Air Temperature

  Var     = "T2M"
  Varname = "2-m Air Temperature"
  Hgt     = 2
  Fx      = 24

  for (Fx in Fx_Hour[2:length(Fx_Hour)]) {
    
    fx_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                             (Height   == Hgt) &
                                             (Fx_Hour  == Fx)  ) %>%
                                      select(-Quarter)
    
    seasonal_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                                   (Height   == Hgt) &
                                                   (Fx_Hour  == Fx)  )
    
    myplot = ggplot(data = seasonal_subset) +
      
      aes(x       = Ens_StDev,
          y       = RMSE_Ens000,
          color   = Quarter) + 
      
      facet_wrap(facets = ~ Quarter) + 
      
      theme_bw() +
      
      theme(strip.background = element_rect(fill=NA),
            ) +
      
      labs(title    = str_c(Fx,
                            "-hr Forecast, CI Seasonal Triangles for ",
                            Varname,
                            sep = ""),
           subtitle = str_c(region_name,
                            sep = "")) +
        
      xlab("Ensemble Standard Deviation") + 
      
      ylab("Root Mean Squared Error") +
      
      geom_point(data    = fx_subset,
                 color   = "grey",
                 alpha   = 0.5) +
      
      scale_colour_manual(values=c("DJF" = "cyan",
                                   "MAM" = "green",
                                   "JJA" = "magenta",
                                   "SON" = "orange"),
                        guide =FALSE) +
    
      geom_point(data    = seasonal_subset,
                 alpha   = 0.7)

      
    print(myplot)
    
  }
## Warning: Removed 332 rows containing missing values (geom_point).
## Warning: Removed 83 rows containing missing values (geom_point).

## Warning: Removed 552 rows containing missing values (geom_point).
## Warning: Removed 138 rows containing missing values (geom_point).

## Warning: Removed 748 rows containing missing values (geom_point).
## Warning: Removed 187 rows containing missing values (geom_point).

## Warning: Removed 904 rows containing missing values (geom_point).
## Warning: Removed 226 rows containing missing values (geom_point).

## Warning: Removed 1128 rows containing missing values (geom_point).
## Warning: Removed 282 rows containing missing values (geom_point).

## Warning: Removed 1312 rows containing missing values (geom_point).
## Warning: Removed 328 rows containing missing values (geom_point).

## Warning: Removed 1484 rows containing missing values (geom_point).
## Warning: Removed 371 rows containing missing values (geom_point).

## Warning: Removed 1620 rows containing missing values (geom_point).
## Warning: Removed 405 rows containing missing values (geom_point).

## Warning: Removed 1844 rows containing missing values (geom_point).
## Warning: Removed 461 rows containing missing values (geom_point).

## Warning: Removed 2008 rows containing missing values (geom_point).
## Warning: Removed 502 rows containing missing values (geom_point).

## Warning: Removed 2196 rows containing missing values (geom_point).
## Warning: Removed 549 rows containing missing values (geom_point).

## Warning: Removed 2332 rows containing missing values (geom_point).
## Warning: Removed 583 rows containing missing values (geom_point).

## Warning: Removed 2528 rows containing missing values (geom_point).
## Warning: Removed 632 rows containing missing values (geom_point).

## Warning: Removed 2696 rows containing missing values (geom_point).
## Warning: Removed 674 rows containing missing values (geom_point).

3b 500-hPa Heights

  Var     = "ISOHGT"
  Varname = "500-hPa Isobaric Heights"
  Hgt     = 500 * 100
  Fx      = 24

  for (Fx in Fx_Hour[2:length(Fx_Hour)]) {
    
    fx_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                             (Height   == Hgt) &
                                             (Fx_Hour  == Fx)  ) %>%
                                      select(-Quarter)
    
    seasonal_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                                   (Height   == Hgt) &
                                                   (Fx_Hour  == Fx)  )
    
    myplot = ggplot(data = seasonal_subset) +
      
      aes(x       = Ens_StDev,
          y       = RMSE_Ens000,
          color   = Quarter) + 
      
      facet_wrap(facets = ~ Quarter) + 
      
      theme_bw() +
      
      theme(strip.background = element_rect(fill=NA),
            ) +
      
      labs(title    = str_c(Fx,
                            "-hr Forecast, CI Seasonal Triangles for ",
                            Varname,
                            sep = ""),
           subtitle = str_c(region_name,
                            sep = "")) +
        
      xlab("Ensemble Standard Deviation") + 
      
      ylab("Root Mean Squared Error") +
      
      geom_point(data    = fx_subset,
                 color   = "grey",
                 alpha   = 0.5) +
      
      scale_colour_manual(values=c("DJF" = "cyan",
                                   "MAM" = "green",
                                   "JJA" = "magenta",
                                   "SON" = "orange"),
                        guide =FALSE) +
    
      geom_point(data    = seasonal_subset,
                 alpha   = 0.7)

      
    print(myplot)
    
  }
## Warning: Removed 340 rows containing missing values (geom_point).
## Warning: Removed 85 rows containing missing values (geom_point).

## Warning: Removed 560 rows containing missing values (geom_point).
## Warning: Removed 140 rows containing missing values (geom_point).

## Warning: Removed 756 rows containing missing values (geom_point).
## Warning: Removed 189 rows containing missing values (geom_point).

## Warning: Removed 912 rows containing missing values (geom_point).
## Warning: Removed 228 rows containing missing values (geom_point).

## Warning: Removed 1128 rows containing missing values (geom_point).
## Warning: Removed 282 rows containing missing values (geom_point).

## Warning: Removed 1312 rows containing missing values (geom_point).
## Warning: Removed 328 rows containing missing values (geom_point).

## Warning: Removed 1492 rows containing missing values (geom_point).
## Warning: Removed 373 rows containing missing values (geom_point).

## Warning: Removed 1628 rows containing missing values (geom_point).
## Warning: Removed 407 rows containing missing values (geom_point).

## Warning: Removed 1848 rows containing missing values (geom_point).
## Warning: Removed 462 rows containing missing values (geom_point).

## Warning: Removed 2016 rows containing missing values (geom_point).
## Warning: Removed 504 rows containing missing values (geom_point).

## Warning: Removed 2204 rows containing missing values (geom_point).
## Warning: Removed 551 rows containing missing values (geom_point).

## Warning: Removed 2340 rows containing missing values (geom_point).
## Warning: Removed 585 rows containing missing values (geom_point).

## Warning: Removed 2532 rows containing missing values (geom_point).
## Warning: Removed 633 rows containing missing values (geom_point).

## Warning: Removed 2704 rows containing missing values (geom_point).
## Warning: Removed 676 rows containing missing values (geom_point).

3c Mean Sea Level Pressure

  Var     = "MSLP"
  Varname = "Mean Sea Level Pressure"
  Hgt     = 0
  Fx      = 24

  for (Fx in Fx_Hour[2:length(Fx_Hour)]) {
    
    fx_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                             (Height   == Hgt) &
                                             (Fx_Hour  == Fx)  ) %>%
                                      select(-Quarter)
    
    seasonal_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                                   (Height   == Hgt) &
                                                   (Fx_Hour  == Fx)  )
    
    myplot = ggplot(data = seasonal_subset) +
      
      aes(x       = Ens_StDev,
          y       = RMSE_Ens000,
          color   = Quarter) + 
      
      facet_wrap(facets = ~ Quarter) + 
      
      theme_bw() +
      
      theme(strip.background = element_rect(fill=NA),
            ) +
      
      labs(title    = str_c(Fx,
                            "-hr Forecast, CI Seasonal Triangles for ",
                            Varname,
                            sep = ""),
           subtitle = str_c(region_name,
                            sep = "")) +
        
      xlab("Ensemble Standard Deviation") + 
      
      ylab("Root Mean Squared Error") +
      
      geom_point(data    = fx_subset,
                 color   = "grey",
                 alpha   = 0.5) +
      
      scale_colour_manual(values=c("DJF" = "cyan",
                                   "MAM" = "green",
                                   "JJA" = "magenta",
                                   "SON" = "orange"),
                        guide =FALSE) +
    
      geom_point(data    = seasonal_subset,
                 alpha   = 0.7)

      
    print(myplot)
    
  }
## Warning: Removed 340 rows containing missing values (geom_point).
## Warning: Removed 85 rows containing missing values (geom_point).

## Warning: Removed 560 rows containing missing values (geom_point).
## Warning: Removed 140 rows containing missing values (geom_point).

## Warning: Removed 756 rows containing missing values (geom_point).
## Warning: Removed 189 rows containing missing values (geom_point).

## Warning: Removed 912 rows containing missing values (geom_point).
## Warning: Removed 228 rows containing missing values (geom_point).

## Warning: Removed 1128 rows containing missing values (geom_point).
## Warning: Removed 282 rows containing missing values (geom_point).

## Warning: Removed 1312 rows containing missing values (geom_point).
## Warning: Removed 328 rows containing missing values (geom_point).

## Warning: Removed 1492 rows containing missing values (geom_point).
## Warning: Removed 373 rows containing missing values (geom_point).

## Warning: Removed 1628 rows containing missing values (geom_point).
## Warning: Removed 407 rows containing missing values (geom_point).

## Warning: Removed 1852 rows containing missing values (geom_point).
## Warning: Removed 463 rows containing missing values (geom_point).

## Warning: Removed 2016 rows containing missing values (geom_point).
## Warning: Removed 504 rows containing missing values (geom_point).

## Warning: Removed 2204 rows containing missing values (geom_point).
## Warning: Removed 551 rows containing missing values (geom_point).

## Warning: Removed 2340 rows containing missing values (geom_point).
## Warning: Removed 585 rows containing missing values (geom_point).

## Warning: Removed 2536 rows containing missing values (geom_point).
## Warning: Removed 634 rows containing missing values (geom_point).

## Warning: Removed 2704 rows containing missing values (geom_point).
## Warning: Removed 676 rows containing missing values (geom_point).

3d 10-m Wind Speeds

  Var     = "M10"
  Varname = "10-m Wind Speed"
  Hgt     = 10
  Fx      = 24

  for (Fx in Fx_Hour[2:length(Fx_Hour)]) {
    
    fx_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                             (Height   == Hgt) &
                                             (Fx_Hour  == Fx)  ) %>%
                                      select(-Quarter)
    
    seasonal_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                                   (Height   == Hgt) &
                                                   (Fx_Hour  == Fx)  )
    
    myplot = ggplot(data = seasonal_subset) +
      
      aes(x       = Ens_StDev,
          y       = RMSE_Ens000,
          color   = Quarter) + 
      
      facet_wrap(facets = ~ Quarter) + 
      
      theme_bw() +
      
      theme(strip.background = element_rect(fill=NA),
            ) +
      
      labs(title    = str_c(Fx,
                            "-hr Forecast, CI Seasonal Triangles for ",
                            Varname,
                            sep = ""),
           subtitle = str_c(region_name,
                            sep = "")) +
        
      xlab("Ensemble Standard Deviation") + 
      
      ylab("Root Mean Squared Error") +
      
      geom_point(data    = fx_subset,
                 color   = "grey",
                 alpha   = 0.5) +
      
      scale_colour_manual(values=c("DJF" = "cyan",
                                   "MAM" = "green",
                                   "JJA" = "magenta",
                                   "SON" = "orange"),
                        guide =FALSE) +
    
      geom_point(data    = seasonal_subset,
                 alpha   = 0.7)

      
    print(myplot)
    
  }
## Warning: Removed 332 rows containing missing values (geom_point).
## Warning: Removed 83 rows containing missing values (geom_point).

## Warning: Removed 552 rows containing missing values (geom_point).
## Warning: Removed 138 rows containing missing values (geom_point).

## Warning: Removed 748 rows containing missing values (geom_point).
## Warning: Removed 187 rows containing missing values (geom_point).

## Warning: Removed 904 rows containing missing values (geom_point).
## Warning: Removed 226 rows containing missing values (geom_point).

## Warning: Removed 1128 rows containing missing values (geom_point).
## Warning: Removed 282 rows containing missing values (geom_point).

## Warning: Removed 1312 rows containing missing values (geom_point).
## Warning: Removed 328 rows containing missing values (geom_point).

## Warning: Removed 1484 rows containing missing values (geom_point).
## Warning: Removed 371 rows containing missing values (geom_point).

## Warning: Removed 1620 rows containing missing values (geom_point).
## Warning: Removed 405 rows containing missing values (geom_point).

## Warning: Removed 1844 rows containing missing values (geom_point).
## Warning: Removed 461 rows containing missing values (geom_point).

## Warning: Removed 2008 rows containing missing values (geom_point).
## Warning: Removed 502 rows containing missing values (geom_point).

## Warning: Removed 2196 rows containing missing values (geom_point).
## Warning: Removed 549 rows containing missing values (geom_point).

## Warning: Removed 2332 rows containing missing values (geom_point).
## Warning: Removed 583 rows containing missing values (geom_point).

## Warning: Removed 2528 rows containing missing values (geom_point).
## Warning: Removed 632 rows containing missing values (geom_point).

## Warning: Removed 2696 rows containing missing values (geom_point).
## Warning: Removed 674 rows containing missing values (geom_point).

3e 10-m Wind Gusts

  Var     = "GUST"
  Varname = "10-m Wind Gusts"
  Hgt     = 10
  Fx      = 24

  for (Fx in Fx_Hour[2:length(Fx_Hour)]) {
    
    fx_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                             (Height   == Hgt) &
                                             (Fx_Hour  == Fx)  ) %>%
                                      select(-Quarter)
    
    seasonal_subset = CI_Ensemble_Stats %>% filter((Variable == Var) &
                                                   (Height   == Hgt) &
                                                   (Fx_Hour  == Fx)  )
    
    myplot = ggplot(data = seasonal_subset) +
      
      aes(x       = Ens_StDev,
          y       = RMSE_Ens000,
          color   = Quarter) + 
      
      facet_wrap(facets = ~ Quarter) + 
      
      theme_bw() +
      
      theme(strip.background = element_rect(fill=NA)) +
      
      labs(title    = str_c(Fx,
                            "-hr Forecast, CI Seasonal Triangles for ",
                            Varname,
                            sep = ""),
           subtitle = str_c(region_name,
                            sep = "")) +
        
      xlab("Ensemble Standard Deviation") + 
      
      ylab("Root Mean Squared Error") + 
      
      geom_point(data    = fx_subset,
                 color   = "grey",
                 alpha   = 0.5) +
      
      scale_colour_manual(values=c("DJF" = "cyan",
                                   "MAM" = "green",
                                   "JJA" = "magenta",
                                   "SON" = "orange"),
                        guide =FALSE) +
    
      geom_point(data    = seasonal_subset,
                 alpha   = 0.7)

      
    print(myplot)
    
  }
## Warning: Removed 304 rows containing missing values (geom_point).
## Warning: Removed 76 rows containing missing values (geom_point).

## Warning: Removed 488 rows containing missing values (geom_point).
## Warning: Removed 122 rows containing missing values (geom_point).

## Warning: Removed 664 rows containing missing values (geom_point).
## Warning: Removed 166 rows containing missing values (geom_point).

## Warning: Removed 792 rows containing missing values (geom_point).
## Warning: Removed 198 rows containing missing values (geom_point).

## Warning: Removed 996 rows containing missing values (geom_point).
## Warning: Removed 249 rows containing missing values (geom_point).

## Warning: Removed 1156 rows containing missing values (geom_point).
## Warning: Removed 289 rows containing missing values (geom_point).

## Warning: Removed 1304 rows containing missing values (geom_point).
## Warning: Removed 326 rows containing missing values (geom_point).

## Warning: Removed 1420 rows containing missing values (geom_point).
## Warning: Removed 355 rows containing missing values (geom_point).

## Warning: Removed 1620 rows containing missing values (geom_point).
## Warning: Removed 405 rows containing missing values (geom_point).

## Warning: Removed 1760 rows containing missing values (geom_point).
## Warning: Removed 440 rows containing missing values (geom_point).

## Warning: Removed 1924 rows containing missing values (geom_point).
## Warning: Removed 481 rows containing missing values (geom_point).

## Warning: Removed 2036 rows containing missing values (geom_point).
## Warning: Removed 509 rows containing missing values (geom_point).

## Warning: Removed 2212 rows containing missing values (geom_point).
## Warning: Removed 553 rows containing missing values (geom_point).

## Warning: Removed 2360 rows containing missing values (geom_point).
## Warning: Removed 590 rows containing missing values (geom_point).